{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXd4VNXWxt89vaRRQkILLXSIlCAdBQQBKQKCIChiodj5LKgI9yo2iuAFEakqRelIlyJgQIoQeqih9wQSSJvJtP39sUmZzJlkUqYkWT+feWTOnLP3mhHX2Wfttd7FOOcgCIIgSj4ybxtAEARBeAZy+ARBEKUEcvgEQRClBHL4BEEQpQRy+ARBEKUEcvgEQRClBHL4BEEQpQRy+ARBEKUEcvgEQRClBIW3DchO+fLlefXq1b1tBkEQRLEiOjr6Huc8OK/zfMrhV69eHYcPH/a2GQRBEMUKxthVV86jkA5BEEQpgRw+QRBEKYEcPkEQRCmBHD5BEEQpgRw+QRBEKYEcPkEQHoNzju8PfI/QqaGQfyFHox8bYdvFbd42q9RADp8gCI8xMWoixu0ch7upd2HjNsTEx6Dvsr7Yc3WPt00rFZDDJwjCI6Rb0jFl3xSkmdPsjqdZ0jB+13gvWVW6IIdPEIRHiE+Lh81mk/zszL0zHramdEIOnyAIjxCsC4ZMJu1y6pev72Fr3EuiIRFvb34boVNDUWVaFYzfOR4Gs8HbZpHDJwjCM6gVanzQ5gPolDq741qFFl90/MJLVhU9JqsJLee3xNwjc3E39S5uJt/E1P1T0W1pN3DOvWobOXyCIDzGhA4T8GXHL1FBXwEMDA3KN8Da59eiQ7UO3jatyFh9ejVup9yGyWrKPGa0GBF9Kxr7b+z3omU+Jp5GEETJhjGGMa3HYEzrMd42xW0cuHkAKaYUh+MWmwXRt6LRpmobL1gloBU+QRBEEVK7bG2HsBUAqOQqVA+q7nmDskEOnyAIoggZGjEUKrnK7picyVFGWwbda3f3klUCcvgEQRBFSJAmCHuG70HT0KZQypRQypRoF9YOe4fvhULm3Sg6xfAJgvBJzsSfwcSoiTh06xDqlKuDz9p/htZVW3vbLJdoVKERjow8gkRDIuQyOQLUAd42CQA5fIIo9lhtVvx28jcsOLoAnHO83ORlvPjYi15fTRaGY3eOof3C9kizpMHGbYhNiMXuK7ux/Lnl6Fmnp7fNc5ky2jLeNsEOt/2NYIxNAdALgAnARQDDOecP3DUfQZRGOOcYtGoQtsRuQao5FQAQfTsaq8+sxobBG8AY87KFBePD7R8ixWyf6ZJmTsPbm9/GM7WfKbbfy9u4M4a/HUAjznkEgPMAPnHjXARRKjl06xA2x27OdPYAkGpOxe4ruxF1NcqLlhWOgzcOSh6/kXwDyaZkD1tTcnCbw+ecb+OcWx69PQCgirvmIojSyu4ru+0KfDJIM6dh5+WdXrCoaAjWBUseV8lVkimPxY341Hh8u/dbDFk9BP878D88ND70yLyeCvK9AmC51AeMsREARgBAWFiYh8whiJJBeV15qOVqWGwWu+MahQbBemmn6YtwzvHj4R/xZdSXiEuNQ7AuGGq5GunW9MxztAotXmv6mtf3Jk7ePYmVp1eCgWFAwwFoVKFRvq6PiYtB24VtkW5Nh9FixB/n/sDXe7/G4dcPo2pgVTdZLWCF0XZgjO0AECrx0TjO+bpH54wDEAmgH89jssjISH748OEC20MQpY2HxoeoOr2qQ5jDT+mHy+9dRnldeS9Zlj8m/zMZn//9uZ10skKmgAwyqBVqmG1mDG40GD/1/Mkhx92TfP7355i0d1LmU5VKrsK49uMwrsM4l8dos6ANDtw4AI4sdyhncvSt3xcrB6wskF2MsWjOeWSe57lTzIcxNgzAKACdOedpeZ1PDp8g8s/+6/vRd3nfTGepVqixasAqPFH9CS9b5hoWmwXlJpdDUnqSw2ctKrXAwj4LUcm/Espqy3rBuizO3juLZnOawWCxV73UKDQ4MeoEapernecYJqsJ2q+0sHFHmWi9Uo+UTx0lGVzBVYfvziydbgDGAnjCFWdPEETBaF21NW7+300cuX0ENm5DZKVIyGVyb5vlMgmGBMl9CACITYjNd8jEXaw7u84hdAYANpsN686twwdtPshzDBmTQc7kkg5frVAXiZ25zu/GsX8A4A9gO2PsGGPsJzfORRClGrlMjhaVW6BllZbFytkDQFltWadhmjrl6njYGsFD40Ocjj+NVFNW9pNcJpdMB2WMubyvoJAp0K9+P6hk9t9Xo9BgeJPhhTPaBdyZpRPOOa/KOW/y6DXKXXMRBFF8UcgU+KTdJ9Ar9XbHtQotvuz0ZeZ7G7dh95XdWBmzEreSb7nFFrPVjFEbRyF0aihazW+F4CnBGL9zPDjn6F+/P+TM8WbKGEO/+v1cnmP2M7PRsEJD+Kn84Kfyg06pQ+sqrTGx48Si/CqSFN9SPIIgSgxj246FXqnHV3u+QlxqHGqXq41pXafhqZpPARChnc6LOiPRkAgGhnRrOt5r9R6+6fwNbiXfwrE7x1AtqFqhwz+f/vUpFh9fDKPVCKPVCACYdmAaKvlXwugWozGt6zSM2TYGDGKlz8Exo/sMhAW6nmFYRlsG0SOiceDGAZy/fx4RIRFoWrFpoex2Fbdu2uYX2rQlCIJzbhc64Zyj/qz6OH//vF1mi06hQ7uwdoi6GpWZyRNRIQKbhmwq0Aav1WZF4LeBdkVsGVQLrIYr710BANxMuon159YDAPrU64NK/pXyPVdR4/VNW4IgiIKQM05+Ov40biTdsHP2AJBmScOOyztg47bM1fiR20cw7I9h2DB4Q77nzciLlyI+LT7zz5UDKmN0i9H5Ht8XIHlkgiB8mmRTstON6JzZLiabCdsubpNM8cwLrUKLaoHVJD9rUalFvsfzRcjhEwTh0zQNzV98W8Zkki0G84Ixhpk9ZtpJN8iYDDqlDlO7TnV5HBu3YWXMSvT+vTf6Le+H9efWe715eQbk8AmC8GnUCjXm9ZoHnVKXmSWjV+oRpAmCgjlGpUP0IajoV7FAc/Wo3QM7XtyBbuHdUCOoBvrV64eDrx1EZKU8w+MAstRLh68bjg3nN2Dt2bV4YfULGLFhRIHsKWpo05YgiGJBTFwMfjr8E24m30TPOj3xRLUn0HpBaySbkmG0GKFgCshlcvSu2xs1gmrghcYv4LHQxzxq495re9FtSTeHjV+dUocDrx5A45DGbpmXNm0JgihRNKzQEDN7zLQ7dubNM5h9eDZ2X9mNe2n3cO7+Oaw6vQqMMfzw7w8Y12EcPm3/qcds3HZxm50eUAYWqwXbL213m8N3FQrpEARRbCmnK4fPOnyGKV2m4ELCBRgtRnBw2LgNaZY0TIyaiEuJlzxmT5AmSLJqWClXIlAd6DE7nEEOnyCIYs+6c+skUyrTLen46bDnVF0GNRoEGXN0q4wx9G/Q32N2OIMcPkEQxR6VXCUpe8DBMePgDFxOvOwROyr5V8LKASvhr/JHgDoAAeoABGmCsH7QegRpgnDlwRVsOr8JF+5f8Ig9OaFNW4Igij2xCbFo9GMju4YpGSiYAiMjR+KHHj94zB6jxYg9V/dALpOjXVg7AMDQNUOx4fwGqOVqmKwmdKjWAWueX1MkHbxc3bSlFT5BEMWe8LLhGBUprc9o4RYcvX3Uo/ZoFBp0qdUFnWp0gkquwhd/f4GN5zfCaDHiYfpDGCwG/H31b7z353setYscPkGUcmzchmWnlqHLoi7o/GtnLDq+CFab1dtm5Zv3W78PtUxaU/5m8k2v9vidEz3HoXGK0WLEouOLJLXx3QU5fIIo5QxdMxSvrX8NOy7vwM4rO/HGpjfw7LJnfaY61FWqBlZF73q9oVVoHT67+vAqnl7yNF5d/yruptz1uG3OKn9NVpNHb67k8AmiFBN9Kxrrzq2zKxRKNadi15VdiLoa5TE7rjy4gqFrhiJkagjqzqyL2YdnF+iGs6TfErz1+FsODUYA0Upx4dGFqDKtCob9MQxmq9nu8xtJNzB2+1h0W9IN43eOx+3k2wX+Pjl5stqTmZLK2WlesTmUcmWRzZMX5PAJohSz68ouB8cHCKf/1+W/PGLDnZQ7aD63OX4/9TviUuNwPuE8Ptj2Ad798918j6WSqzC5y+RcC5ws3IKVMSsxYdeEzGPH7xxHg1kN8P3B77H14lZM2TcF9WfVx9l7Zwv0nXLyfbfvEaAOgFouQk5KmRJ+Kj/M7jm7SMZ3FXL4BOEFbiXfwrJTy7A1dqtkn1RP4ay9oEahQXldeY/Y8P2B75FiSrGLZaeZ0zDvyDzEpcYVaEyp1XR2DBYDfjz8Y+b7UZtGIdmUnNlbN92ajqT0pCLbVK1bvi5Ov3kaY1qNQafqnfDW42/h5OiTLmv0FBUkrUAQHmb8zvGYun8qlDLxKK9VarHjxR1eKbvvX7+/pFOTMzkGNRrkERv2XN0j2cRcLVfj5N2T6Fyzc77HfKL6Ezh8O/cU7+T0ZHAuqnL/vfGvw+ccHLuu7LI7djHhIib/MxmHbh1CowqN8HG7j9EguIFLNlXyr4RvnvrG9S/hBmiFTxAeZNvFbZh+YDqMFiOSTclINiUjLjUOPZb28Gi2RgaBmkBsGbIFwbpg+Kv84a/yR1ltWawbtA4V9BU8YkPtcrUli6ZMVhOqBUnr0+fFmy3ehCwP99asYjMwxiBjMqgU0k3Us+fIH79zHE3mNMHCowtx9M5R/HbyNzw+73Hsvba3QDZ6A3L4BOFBfjr8k2QLvYfpD/HvTcdVpidoG9YWt9+/ja1Dt2LLkC248/6dAq2qC8r7rd+HWmGfTqmWq9G6SmuElw13OP/qg6vYGrsVVx5ccTpmjTI1MKb1GMkbiQwy6JX6zEIsxhheeuwlaBQau/O0Ci1eb/Z65vv3tr6HFFMKLFyE4KzcilRzKt7Y9IbL39XbUEiHIDyIs05MjDGkmhxvBJ5CLpOjddXWXpm7cUhjrBm4BiM2jkBcahw45+hVpxcW9Flgd57JasKQ1UOw8cJGqOVqpFvT0T28O37v/7vDDQMApnSZgqdqPoVp+6fh8oPLUMqUUMqUaFqxKT5u9zHqla+Xee70p6fjUsIl/HP9HyjlSpitZnSp2QUTO07MPGf/9f2S9p+KOwWz1ezRbJuCQtIKBOFB5kbPxZitYxwkdHVKHeI/jC+SMvuCwjmHyWqCSq5y6Cvrqfnvpt6Fn8oPfio/h8/Hbh+Lmf/OtCtg0iq0GBk5EtOfnl4kNpyJP4Pz98+jQXAD1C5X2+6z0KmhuJvqmMOvU+qQ8kmKV36zDEhagSB8kGGPDUNESAT0Sj0AsTmqU+jwU8+fPOLsOeeYfXg2wmeEo+yksui7rC/Oxp/FpL2TUG5yOei+1qHa99WwPGZ5kc9ttBgx+9BsdPy1I55d9iy2xm61+5wxhlC/UElnD4ibZc5qVYPFgPlH5heZjfWD66NPvT4Ozh4A3mn5jsN/o4ywjzedfX6gFT5BeBiz1YzVZ1Zj3VmxMfp689fRqEIjj8z9wbYPMPvw7MwnDAYGpVwJOZPbOVOdQocVA1bgmTrPFMm8JqsJbRe2xen405lz65V6/F/r/8MXHb9waQzlRKVkCquMyWAZb3G707XarHhj0xtYdGKRCClZ0tGnXh/8+uyvkiElV+GcI9WcCo1CA4WsYFF2V1f45PAJopSQaEhEpWmVJHXjpWga2hRHRh7J9RzOOVafWY3Zh2cj1ZSKwY0GY0TzEdAq7eUNFh9fjNGbRjtsWKvlalx+9zIq+ufeg/ZG0g20XtAaN5JuOHzWukpr7Ht1n0vfqSiIT43H+fvnUbNMzTztzoutsVvxxuY3cO3hNShlSoxoPgKTu0yWrI3IDWpxSBCEHWfvnYVarnbZ4bvSKeqdP9/Bz0d/znTkJ+6ewJKTS7DvlX12m5gbzm+QzE5SyVXYc20PBjYc6HSOmLgYtFnYxmHfQyFTQKPQYFaPWS59n6IiWB+MYH1wocc5dPMQ+q3ol/m9LDYL5kbPRaIxEb8++2uhx5eCYvgEUUoICwyT1It3Rl4FRZcTL2P+kfl2jtxgMeBs/FmsPrPa7twQfYhkiiQgqn1z450/30FyerJDOKe8tjxOjj6JphWb5nq9r/LVnq9gMDvuSayIWYH7affdMqfbHT5j7APGGGeMeaZOmyAISSoHVEa3Wt0c8s2VMqXDMZ1Sh286514VGnU1CgrmGCRIMadgS+wWu2Mjmo+QDFPoVXp0rN4x13n2XN0DDsfQc1xqHMICw3K9tqDcT7uPc/fOSeoMFRVn752V/F4quQrXHl5zy5xudfiMsaoAugBwj/UEQeSLpf2XYlDDQVDL1VDJVagWWA3rB6/HnJ5zULNMTWgUGjQLbYb1g9bjiepP5DpWsD5Ysn+rUqZEJb9KdscahzTGnF5zoFfqEaAOgJ/KD2GBYfjrpb8gl0mv/DPQq/SSxzVKTZ6aOfklOT0Z/Zb3Q+VplRE5LxIVplTAz0d/dulazjlO3D2Bo7ePulQ13aJyC6cVxrXK1sq37a7g1k1bxtgqABMBrAMQyTm/l9v5tGlLEJ7BaDEi1ZSKstqyBc5uMVvNqDq9qiiWyrZS1Sl0ODH6hKTTSjOn4eCNg/BT+SGyUqRLc0vl32sUGrzW9DXM7DGzQLY7o+dvPbH94naYbFnaPjqlDusHrc+1+jj6VjT6Lu+LRGMiAMBP6YcVA1agfbX2Tq85d+8cIudGIsWcpZWvU+rwVou3MKnLpHzZ7fU8fMZYbwA3OefH3TUHQRAFQ6PQoJyuXKFSGY/eOYrONTrDT+UHjVwDf6U/AtWB+P25352uUHVKHTrW6IgWlVu4PPfEThPRPbw7NAoNAtWB0Cq06FyjMyZ3mVxg23Ni4zaM+XMMNl3YZOfsAXGTmvSPcwecYkpB50WdcT3pOlJMKUgxpeBO6h30WNoD99Kcr3Hrlq+LPa/sQecanaFX6hEWGIZJT03Ct099W2TfKyeFytJhjO0AECrx0TgAnwLo6sIYIwCMAICwMPfE4wiCKFom7JqA7/Z/B6PFCM45NAoNOtfqjOXPLc93SmFeqOQqrH5+NS4nXsaZe2dQp1wdSY2dwjDx74n4Kfonp5/nFlNfc2YNrNyxa5WVW/H7yd/xdsu3nV7bJLQJdry0I3/GFoJCOXzO+VNSxxljjQHUAHD80V28CoAjjLHHOed3cowxF8BcQIR0CmMPQRCFJzYhFgmGBESERDhs5gJCInjKvil26Z0GiwHbLm7DsTvH8Hjlx/Oc49idYzh59yTCy4ajVZVWLq32a5SpgRplauTvy7iAjdsyb15SKGSKXDeW41LjkG5xzH4yWAy4k3JH4grv4ZY8fM75SQCZ2qqMsStwIYZPEIT3uJV8C32W9UFMXAyUciVs3Ibvn/4erzZ71e68zRc2QyK5BAazAevPrc/V4RstRvT6rRf23dgHGZOBc4465ergr5f+QhltmXzbzDnHz8d+xjd7v0F8ajxaVWmFyV0mIyIkwuUx0sxpDpIN2fFT+eGT9p84/bxDtQ5CcM1mn9Hjp/LDk9WfdNkOT0B5+AThw1hsFiQaEj2ild9jaQ8cvX0UBosBSelJSDGl4J0/38G+6/ZVrBqFRjKzRiFT5KkHNGHXBOy9vhdp5jSkmFKQak5FTFwMRm0aVSCbv4j6Am9veRuxCbF4mP4QWy9uRduFbXHu3jmXx9Ar9U61/4N1wTg68miu6Z8tKrVAt1rdMvWRALFX0bJyS4/KTLuCRxw+57w6re4JwnVs3IYJuyagzKQyCP0uFKFTQ/HL0V/cNl9MXAwuJFxwiEUbzAZM3z/d7obTt35fyRuQQqZw2iXrUuIldF3c1SEUBAAmmwlrz6yF2WoG5xzzjsxD7Rm1EfRtELov6Y6Td09KjplqSsXkfyY7VOAazAZMjJooeY0UjDFM7TrV4WalU+iweuBqVA+qnuf1KwaswMzuM9G6Sms8XvlxTO0yFZuHbJZMW/UmJK1AED7If3b9B9MOTMt0ZvFp8Xhzy5sI0gbh2XrPFvl8calxmS0Xs8PBsebsGii+UKBZxWaY1WMWWlZpid/6/4Yha4ZAzuTg4LDYLJjVYxZqlqmZeW1SehKWnFiCI7eOYHnMcklphQxs3AYrt+LL3V9i6v6pmd/7z4t/Yu/1vTj8+mHULV/X7ppLiZck89it3IqDNw/m6/sPbjQYgepATNg1AZcfXEbjCo3xdeev0aZqG5eul8vkGN50OIY3HZ6veT0NiacRhI9htppRdnJZpJhSHD5rXKExTow+UajxN57fiCn/TMHtlNvoWqsrPm3/KfRKPUK/C81TZ0ev1OPIyCOoU64OktKTsOXCFlhsFnQL74ZyunKZ511MuIhWC1ohzZzmsAKXomXlltjx0g5UmFLBIZ4uZ3IMbjwYi/sutjueYEhApe8qScpFdAvvhi1DtjgcL6l4PQ+fIIiCkZSe5LSkv7Al99P2T8Pzq55H1LUoXEi4gLnRc/HYT48hzZyGL578wi4OLYXRYsTUfVMBAAHqADzf6HkMiRhi5+wBYOTGkUgwJOTp7JUyJQLUAZjXax4uJV6S7Bpl5Va7JuNxqXF4df2rqD2zNmRM5iDvoFPq8Fn7z3Kd917aPVx9cBW+tOD1BBTSIQgfo4y2DPzV/khPc1y5Ng5pXOBxU02pGLdznN0q3mwz46HxIabum4rvnv4OjUMa4/sD3+Ni4kVcf3jdYfVs5VYcu3Ms13msNit2X9md50azDDJ0qdUFC3ovQKhfKO6n3ZdMbwSQ2ZAkzZyGFvNa4Hby7cysGDmTQwYZFHIFKugr4IfuP6BtWFvJceJS4/DC6hew99peyJgMZbVl8XOfn9GlVpdcbS0p0AqfIHwMGZPh685fS24iTnoqfyX32dlxaYdkyMZsM2P7pe0ARCjkz6F/IurlKMnceIVMgeYVm+c6D2PMpc1KG2yIuhqFlvNa4sL9CyinK4eBDQdCq7DX0tcpdRjXfhwAYOmJpbifdt8uBdLKrdAoNdj8wmZce+8a+tTrIzkf5xxdF3dF1NUopFvTYbAYcDP5Jp5d/izO3z8vec21h9dw4u4Jt4qoeRJy+AThQXZf2Y32P7dHyNQQdP61s0PKYwavN3sdi55dhIbBDRGgDkDbqm2x7cVtLm8iSrHq9Cqnn4X4hdi9r+hfEQMaDIBOYX/T0Sg0+KDNB7nOI2My9KvfT3ITOCcpphRcT7qOZ357RmTo9JqH4U2GQ6PQQCVToWpAVSzrvyyzwfqBmwecbv5efnDZ4SZltpoRlxoHq82KI7ePIDYh1iFf3mQxYea/9po8d1LuoM2CNqj7Q120W9gOFaZWwIqYFXl+H1+HQjoE4SE2nd+EgSsHIs0i4to7U3fiwOID2DB4AzrV6ORwfv8G/dG/Qf8imz8mPsbpZ73r9nY4tqD3AlQPqo5Zh2YhKT0Jrau0xozuM1xScpzVYxZOxZ3C1YdXJTefs8PBcSv5Fk7GnURESARmPTML07tNR6opFUGaIDsnXq9cPWgVWsmN3RpBogr3bspdrDq9CuvOrcM/1/+BjdugVWjxXIPnJOsHLNyCiwkXAYhw1L83/8VLa1/C5cTLsMIKI8RT0fB1wxFeNhzNKjbL8/v7KpSlQxAeovbM2ohNiHU43iS0CY6OPOr2+YeuGYrfT/3uEFtXypS4/O5lVA6oXKTz2bgNu6/sxrzoeVh7dm2uzVcC1AHY9MImtAtrl+uY99LuIXxGOB6mP8w8ppApULNMTZx58wxWn16NYX8Mg9lqhoXbN0zRKrROj3/R8Qu0rtIafZf3Rao5VXKzWcZkGBox1G3dqAoDZekQhA9htVklnT0AnIo75REbxrYd66CNo1Vo0btu7yJ39oBwkJ1qdMLvz/2ODYM3oGvNrgjVh0o2TeGcI7JSnv4K5XXlETU8Cs1Cm0EpU0IpU6JLzS6IejkKD40PMeyPYTBYDA5OHRDaNjqVzm5vRClTooy2DAY0HIDuS7sjPi3eaWaRjdtw7UHxbu1BDp8gPIBcJkcZjbRWTIg+RPJ4UdM4pDE2Dt6IuuXqZvaDHfbYMCzpt8Ttc3ep1QVbX9yKi+9eRP3g+pnpn3Imh06pw4/P/Cgp1CZFREgEokdG495H9/Dg4wfYPGQzQvxCsOnCpjybqVhsFkx/ejoaBjdE1YCqGNl8JI6MOIIdF3fkmVWkVWjRLbyba1/YR6EYPkF4iI/afoSJURPtVpB6pR6ftv/UYzZ0rNERZ986i1RTKtQKNRQyz7oAnVKHg68dxNKTS7Hx/EaE+oVidORoPBb6WL7HClAH2L232hwlinPSqEIjjGg+AiOaj7A7ft9wP9eaAQaG8rryGBk5Mt92+hLk8AnCQ3zU9iOkmlIx/cB0cHDImAwft/0YoyNHe9wWZ20DPYFWqcVrzV7Da81eK9Jxe9TuAcsmx1BO5rwKLaZ0mSL5WacanST7y2agkqtwdORRBGmCCm2nN6FNW4LwMOmWdMSlxiHEL6TIm4WUdhYeXYg3N78Jq80Ki004f5VchccrP45vOn/jtCALAAK+CUCyKVnys7DAMFx976pbbC4KXN20pRU+QXgYtUKNqoFVvW1GieSVpq+gY/WOWBGzAgaLAX3q9kHTik3zvO5m0k2h53PZsfsUA8PYtmPdYa7HoRU+QRA+A+cce67twe4ruxGsC8bzjZ5HWW3ZAo1lMBtw6NYh+Kv80SS0idOuWoduHkKnRZ1gspgc+tkCwIsRL+LXZ38tVP9fd0MrfIIg8sRqsyLdmp5n4xJPYLFZ8OyyZ7H7ym6kmdOgVWjx0Y6PsPmFzWhfrX2+xlpyYglGbxoNGZPBxm0I0Ydg85DNqFOujsO5w9cNdygOY2CICInAhsEbStTTGKVlEkQpxGQ14d0t78L/G38EfBOAOjPrYMclzzXTlmLx8cXYdWUXUs2p4OBIs4iuWM+tfM6lDJwMTtw9gZEbRyLFlJLZuetS4iU8tegph9TLREOipI4OB8fNpJslytkD5PAJolTy6rpXMe/IPBgsBli5FRcSLqDPsj44cvuI12xaeGyhZGqkwWzIl12zD812UN3k4HhgfICoq1F2x3PbNNcqtU4f8HZHAAAgAElEQVQ/K66QwyeIUkZ8ajxWnV7loEdjMBvwzZ5vPGLDxYSLeGvzW+jwcwd8uO1D3Ey6CQbnMfKM+DnnHJcTL+NOyh2n595JuePQqjGD+2n37d7rVXo8Hf60g9CbVqHFqMiC9dn1ZcjhE0Qp4+rDq1ApHFe2HByn7512+/wHbxzEYz89hjnRc7Dn2h7M+HcGGv7YEN3Du0s2YNEr9WhWsRkWHV8Ev2/8UGtGLVT6rhIazmqI6w+vO5zfq24vyXHMVrOkVs/C3gtRP7g+/FR+8FP5QavQ4unwp/Fhmw+L5gv7EOTwCaKUEV42HCarYzaKnMnxeKXHCzwu5xw//PsDanxfA4HfBuKZpc/gdLzjDWTkxpFINadm5smbrCYkpSch6moUutTqAr1SDzmTQ6/Uw1/ljzXPr8HOSzsx7I9hSDOngT/65/S904icF+kQl3+h8QsILxtup6uvV+rxUduPHGSgASBYH4xjI49h29BtmNtzLqJHRGPt82slu28VdygtkyBKIe9vfR8/Rf9kFzP3U/khekS0ZCaLK/zf1v/DnOg5mWMyMPip/HB81HHUKCOki9Mt6dB9rZPUrVHJVLjxfzcQmxAr0jL1wRjQYAACNYGoM6MOLiRecLiGgWHHSzsc5KXTzGlYcGQBlscsR5AmCG+2eBPda3cv0PcqDpBaJkEQTpnSdQq+6vQVqgRUgU6pQ+canbF3+N4CO/tEQyJmH55tdwPh4DCYDfj2n28zjylkCqcbpSabCWHfh+HC/Qv4pP0neK3ZawjUBAIAriVJq1RycJyOc3yKuJd2Dzsu78DhW4ex+8pu/HH2DySnS1fRliYoD58gSiEyJsN7rd7De63eK5Lxzt0/B7Vc7dBC0cItOHjjYOZ7uUyOFyNexOITiyXbLRotRozaNAoRIRE4FX8KGoUG3cO7o5y2HG6l3JKcu0O1Dnbvk9OT0WJeC9xLuwcbtyHdmo5fj/+Ko3eO4uBrB326gMrdkMMnCKLQVAusJunAZUyGeuXr2R37vtv3uJl8E9svbndoNwiImH6L+S2gUWjAwMDBMbr5aEzdP9VB4KxWmVqICI2wO7b4xGKkmFLswkbp1nScuXcG/1z7B8mmZCw8thBmqxlDI4aiX/1+LvXgLQmUjm9JEIRbqehfEb3q9HJoQK5RaBx0aHRKHTa9sAn/1+r/IGeO+vVWLoTPUkwpSDYlI8WUgh8P/4hP230KpUwJ9uiflpVb4sSoEw7XH71zVDKf38ZtGL97PAasHJDZAvHlP17Gcyuegy/tZboTcvgEQRQJi/stxosRL0Kj0EApU6JmUE2sfX6tU/Gy3vV655p7nx0Zk6FRSCMkfZKE02+eRsLYBBx47QB0KkdJiIgKEU6lIvZd32fXBD3VnIptF7fh76t/u2RHcYccPkEQRYJGocGcXnPw8OOHiP8wHrHvxKJrra6S5564ewLdlzpmzUit+AGxOjeYDdAoNKhXvl6uuvQvPfYSdEqdXZhGJVeJhugSN5g0cxo2X9ic19crEZDDJwiiSFHJVQjUBOa6OfrKuleQlJ5k13uWgTl15BabBU+HP42k9CSsjFmJ5aeW44HxgeS5gZpAHHj1ADrX6Aw5k0MlV2Fgw4H4rP1nDhW1AKCUK522n8yuq18ScGsePmPsbQBvAbAA2MQ5/yi38ykPnyCKDw+ND7E8ZjluJt1Em6pt0KVWF5c2P1NMKSgzqYzLjlTGZPi287eoU64OXljzQuZTgMVmwbze8zCk8RCn12b4N8YYktKTUHlaZQdlTK1Ci7NvnUVYYFjmsfjUeIzeNBrrz62HjdvwZPUnMbfXXNQsU9Mlmz2N1+WRGWMdAfQBEME5T2eMVXDXXARBeJYjt4+g468dYbFZkGZOg5/KDxEhEfjrpb/ybEaukClcjt1nnN8tvBtazm/poP/z+vrX0T6svZ2zzk72p4wAdQA2vbAJfZb1yczgsdqsWNR3kd31VpsV7X5uh0uJlzJvSruu7EKr+a1w6d1L8FP5uWy7r+HOkM5oAN9yztMBgHMe58a5CILwEJxzDFw5EEnpSZnZMCmmFBy9fRQzDs7I83qNQoNu4d1cbqCulCmx9uxayRCRlVux/NRyl23vUK0D7n5wF2sGrsHy55Yj7sM49Kvfz+6c7Ze241byLbsnEBu3Ic2cht9P/u7yXL6IOx1+HQDtGWMHGWN/M8ZaSJ3EGBvBGDvMGDscHx/vRnMIgigKLj+4jFvJjkVQBosBvxz7xaUxFvRegDrl6sBf5Q+9Ug+dUodqgdWgkjlW4QZpgqBX6iU18S1Wi0OIJi9UchU61+yMbuHdJLN5zt8/D7PVsT4g1ZyKmPiYfM3laxQqpMMY2wEgVOKjcY/GLgOgFYAWAFYwxmryHJsGnPO5AOYCIoZfGHsIgig4nHPcN9yHVqGFXuWoNpmBKzLGeRGsD8ap0afw99W/cTnxMh4LfQw1gmogcl4k7qTcQZo5DUqZEkq5Er88+wuqBlTF+F3jHcZRypUI9QtFuiUdaoXapbnzolGFRlDKlUi32mvq+yn90DQ07/64vozbNm0ZY39ChHR2P3p/EUArzrnTZTxt2hKEd4i6GoVX1r2C60lCbrhXnV6Y33u+ZNYM5xz1Z9XHufvn7I7rlDp88eQXeL/N+wW2I9WUisUnFmPn5Z2oHlQdoyJHZW6Ujt0+Fj8c+gEGsyGz4lYpU2Y2Kln+3HJ0C+9W4Lkz4Jyj+dzmiImPyVQVVTAFQv1Dcf6t8z7ZGMXVTVt3OvxRACpxzicwxuoA+AtAWM4VfnbI4ROE57lw/wKazGliV52qkqvQolIL7H1lr+Q1J++exBO/PAGzzQyjxQi1XI1WVVph85DNuXaRKiz/XPsH84/Mx5KTSxyyfHRKHS69c0lSAjm/JKUnYez2sfjt1G+w2CzoU7cPvuv6HSr6Vyz02O7A61k6ABYCWMgYOwXABGBYbs6eIAjv8L+D/3PQxzdZTTh65yhi4mLQsEJDh2sahzTG9THXsebMGtxKvoU2VdugXVg7twuTtQ1rizP3zmDl6ZUODt/GbVh2ahnebfVuoecJUAdgds/ZmN1zdqHH8iXc5vA55yYAQ901PkEQRcPZe2clc+KVMiWuPLgi6fAB0R7wxcdedLd5DiQaEiUbuKRb0pFgSHB6Hecc2y5uw5qza+Cn8sOwx4YhIiTC6fklEaq0JYhSTvuw9tDIHXPn063pPukQu9TqItmNSqfUOZVysHEbBq4ciP4r+mNu9Fz878D/0Gp+K8z6d5a7zfUpyOETRCnnjRZvwE/tZ6djo1PoMKjhIFQNrOpFy6RpEtoEAxoMsOtbq1fq0T28O9pUbSN5zdbYrdgSuyVTOM3KrTBYDPhg+weITy096eDk8AmilBOsD0b0iGgMajQI5bXlUSOoBr7s9CXm957vbdOc8nOfn7G472L0rNMTPcJ7YGGfhVg+YLnTPYSVp1faqWRmoJQpsf3Sdneb6zNQAxSCIBAWGIYl/ZZ42wyXYYyhb/2+6Fu/r0vnZ6hnSvXSzanhX5KhFT5BECWe4U2GO9X4eTr8aQ9b4z3I4RMEUeJpXqk5Pn/yc2jkGvip/OCv8oefyg/rBq1z2iylJOJWeeT8QoVXJYzYWGDrVsDfH+jTBwgM9LZFRCnndvJtbLu4DXqV2OTNTUKiOOH1StuCQA6/BPHRR8APPwCcAwqF+Pe6dUDnzt62jCBKHK46fArpFDfu3QMuXgSsjsqBPsPu3cCPPwIGA2A0AikpQGoq0LeveE8QhFcgh19cSEgAuncHqlQBIiKAihWBtWu9bZU0P/8sHHxOGAP++svz9hAEAYAcfvGhTx9g504gPR1ISwPi44GhQ4EjR7xtmSMmx7L3TMyOOuMEQXgGcvjFgfPngehoR0dqNALffecdm3LjhRcAvcRmmMVCMXyC8CLk8IsDN28CKgnJWZsNuHzZ8/bkxTPPAD17CqfPmLBdqwXmzRMZOwRBeAWqtC0ORESIUE5O1GrfXDHLZMDvvwN79gAbNgABASL8VKOGty0jiFINOXxPkJICLFkC7N8PNGgAvPIKEBzs+vXlygHvvQfMmCHi94BIdQwIAN55xz02FxbGgA4dxIsgCJ+AHL67uX0baNECePBAZK5oNMDXX4vVb0Q+pGe//hpo1AiYNg24f19k7Iwf7/qNIy0NWLMGuHEDaNkSePJJ4ZQJgig1kMN3Nx9/DNy9KzYsAbHRajQCw4eLjVhXYQwYMkS88suZM2KlbTSK3HitFmjSBNi+XdwItm0DlErg6acBP7/8j08QRLGAHL672bAhy9ln58QJEerxhIMdNEg8FWRUVaekAIcPi2yaLVuEswfEJvCqVUC3wjeCJgjC9yCH726ksmsAsWJXeODnv3VLpHXmlNAwGoE//hDHs1e/9u8vsoKCgtxvG0EQHoXSMt3Nq6+KuH12lEqxis553B3kppUk9ZlMJjRvCIIocZDDdzfjxwPt2gE6nchL9/MDatcGFizwzPyVKwM1azoel8sdjwFCoycjE4ggiBIFhXTcjUYjNkePHAGOHxfOt0MHz2bILFsm5jSbRaaQn5/Q5Ll2zdG5cy4ygPKLwSCygK5cASIjgS5dxNMCQRA+Azl8T9GsmXh5g8aNhXNfsQK4fl2kZXbtKjKF1qwRNwGZTNycPvoIqF49f+PHxgJt24qbR1qaeJqpWxf4+29piQWCILwC6eGXZjgHduwQTwBqNfDSS0CrVvkfp1Ur4N9/7fcE1GpRLPbtt1nHTCaRBfTXX0BYmNjfqFKl8N+DIEo51ACF8AyJiUBoqLRCZuXKotALEKmgbdoI7Z+UFHFDUCiAjRtFERhBEAWGGqAUd4xGYNYsoH17IUa2cWPuGTf5ITZWaOmfOFG4cXbvBkaNkq4zAOztnTYNuHBBOHtAaAOlpopCMputcHYQBOESFMP3RUwm4ehPn87aVN25Exg5Evj++4KPazaLYquNG0V9gMUCNG0KbNqU/36zn34qtH2kGp0AYgU/dGjW+2XLpLtdPXwInDsH1K+fv/kJgsg3tML3RVauFHII2TNojEbgf/8D/vOfgo/7zTfCuRuNQFKSGP/QIeCNN/I3zqVLwPTpzp29n59w4OPHZ7Vi1Gqlz7XZnH9GEESR4jaHzxhrwhg7wBg7xhg7zBh73F1zlTg2bXLuTL/9VqR5FoSffhLpk9nJ2EjNrUtVTrZvd55y2bKlkEYePBioVk3E6WvUEMd1OvtzZTJRk5DfrCCCIAqEO1f4kwF8zjlvAmDCo/eEK1So4LwwymQSq+v8YjIJxU4pbLb8tR7095d2+EqlSPc8dgz4/HPRhxcQufm//CKyebRakarp7y82e9esye83IQiigLjT4XMAAY/+HAjglhvnKlmMHJm7zs7du/kb7+RJkf7ozKnXquVc80eKXr2kjysUYmU/ebJjQZfBIDZso6NFaGrZMuDqVTE3QRAewZ0O/z0AUxhj1wFMBfCJG+fyPcxmsXr98sv8h0zq1xftAKVQq0X7QFfhXJwfH++YTZNR7XvtmniqmD/ftTH9/cXGb2CgaMISECBW7l9/DcTFOb+xXLggvturrwI9ekjf1DgXm8GVKoknhmbNRAEXQRCFh3Ne4BeAHQBOSbz6AJgBoP+j8wYC2OFkjBEADgM4HBYWxksEd+9yXrMm5/7+nDMm/l29Oue3b+dvnJkzOddoOBduUPw5LIzzhATXxzhyhHM/v6wxcnvpdJxv3er62EYj51u2cD5/PudNmnCu1XKu14vvLDV+mzZ5j/nf/wo7ctp14IDrdhFEKQPAYe6Kz3blpIK8ADxEVmEXA5CU1zXNmzd33y/iSQYN4lyhsHdaCgXn/fvnf6yoKM779uW8VSvOJ07Mn7PnnPN9+zgPCHDN4QOcP/lk/sa32TgPD+dcLs/7ZrJ7d+5jGQzihiF1fdeu+bOLIEoRrjp8d4Z0bgF44tGfOwG44Ma5fIu1ax3DJxYLsH59/oun2rcXoaH9+4HPPgPKlJE+z2YTCpxNmoi4+Pvvi6YnzZvnT6jt2rX82bdvH3DnTlb6ZQaMidCPUilaOf7xB/DEE9JjZHD7tvPPTp3Kn125ERMj+gq3bg18+KHoGUAQpQB3Fl69DuB/jDEFACNE6IZwF2+8IRqlZ6Rz/vADsHq12LBdvBh4/nlx0zGbRZaMweBY4SqX5+2Uc3LnjnTGDudAp07C0ecko9ViUJD9zSg01HnVbb16+bPLGX/9BfTuLSp9rVahYjp/vugARhvIRAnHbSt8zvleznlzzvljnPOWnPN8NHAt5jz7rOOGpEIhslsKK4ucnp4lT5DB9evAr7/a5+6bTGKj9uefxbynTgklzFdfFSmS06bZ58XLZKJg6rPP8mdPy5bCppzodKJHbnaSk0Wlb2AgEBKSpaiZgVYLvPOOY76+TifSPAsL58CIESKDKOOJxGQSRWgff1z48QnC13El7uOpV4mJ4d+9y3mNGlmbtn5+nFerlv9N2+zcvy/2AFQqsR/QtKnYkOWc89Wrncfpe/Z0Pua6dWJvICyM8xdf5Dw2tmC2vfWWfexdrRZx/ZQU+/M6dRKf5Yztnz2bdY7VyvlXX3Fepoz4vG5dsTFcFNy/L34/qd+pTJmimSMvzp7l/OWXOW/UiPPnn+f8+HHPzEuUaODtTduCvEqMw+ecc5OJ85UrOf/iC85XrOA8Pb3gY9lswsErlfZOyt+f81u3RAaLVCaOQsH5u+/mby6rlfNZszivV4/zSpU4f/11MUde9i1dynnLlpzXr8/5+PGcJyban3PunKOzz7Bx1CjnthQlBoO0DYDIqnI30dHiv1PGBrdM5tpmNkHkgasOn8TT3IVSCTz3XNGM9e+/ohF5zvx2kwmYOxeYMEHoy58/b79ZrFLlrZNjMolN5mPHRIglKgpYvjyrcOrnn8Vm8+nTQNmy4hjnwIYNQs3zyhVR1DV4MLBrl3NdnFmzpEM/FovQDZKiqDtmaTSiSfvq1fa26HTAu+8W7VxSjBljH46z2cTv/MYbYiOZINyNK3cFT72K1Qo/IYHzzZs5P3hQrHBd4cQJsXpeuVKsNl1l6VKxmpdamQ4cKM65dYvzyEixasxYPb7wQu5PFnFxYmWb8XSg1UrPodVy/s03WdeNGeOYKy+TiSeCmzel53G2slYqOf/kE9d/i8KSnCxSPLVazgMDRW3DG28U/dOEFM5+A8bEEyFBFBBQSMdNWCwi3i2TiXiwTsd5rVqcX7zo/BqrVThfnU44GH9/zsuV4/zkSdfmPHlS2hnrdJx/95045949zsuXz3L4GY66b1/n4w4b5hgmcvbKyIO/csW+GCyn05eqNViyxHnxl0rl+t7GjRsiXPT885z/8APnSUmuXSfFxYuc79wp9ls8RUiI9G+g17u+aCAICcjhuwOzmfPmzaX/p61b1/n/tL/84rgizogbu/o/+jPP2Dt9uZzzChVErHzLFnHTkapw1WpF/FwKZ08NOV9yOedvvimu+fXX3Ct31WrHeZYvdz7XsGGuff+MfYqMVbJOx3mVKp512IXl228d/x5otZy//763LSOKOa46fNLDzw+rVolYtxTXrzvvIDVnjqOYGCBE0E6fdm3uNWtEkVBoqEhrfP55IUT2yy8iLn3xonAhOVEqndvlaozcahX59AcPAuXK5X6d1Gfdu0vn1+t0wFtvuWbDsGEi/p0Re09LE7/fhAmuXe8LfPihKPjSaMR/Q40GGDBA9CkgCA9ADj8/LF/uWFGaAefO5YedCacx5rossUolctFv3xbzLF0qqm7HjZO+mWRgsQDh4dKfDRrkukrmzZvAU0+JSl61WvocpVLciHLi7y+auuh0ItdfpxPObtw4IDLPNpyinuDyZcfjZrN0YZevIpMBM2eK33LbNlHV/Ouv4ncjCA9QMrJ0rl8H9uwBgoNFdaczLfnC4ufn/DPOgRYtpD8bOlSs5HM2H9HpgMaNReVpdLRwjI0bu16cdfJk7jLKarVw0E2aSH8+aZL43Vx9yrBagRUrRLVq9+5CkiDjqUKtFjeWadOkr804f/168Tt06yYyi1zB2Q0GKJ7dssqWBR6nfkCEF3Al7uOpV75j+DabyBjJ2Aj19xeZItkLeYqS3bulY/EA57NnO7/OYBAFThmxb41GjPPXX2JD099fFE7p9ZzXrs35+fOu2XPpkvPMGsY4HzyY8wcPch9jwwbn30nqNWaMuM5qFXH1yZPFa+vWrEwXo1Fk6xRl5knXro6CdDod55MmFd0cBFFMQanYtF2zxlFdkTHhNN2V9fDf/4qNQ7VabGYqlUIeOC8sFs7/+IPzt98WKY43b4oqy5zOljFR+epqmmC7do6ZNjod59u3u3b9vn2uyyf7+Ymq3vv3Of/oI1FN27y52Mi12YTN48aJ/yYZaY/TphXNf4s7dzhv0EDY4OeXlYFE6YwEUUocfufO0o5Jrxc57+7ixg3h5NauzV8+fU5GjZKWFc7IoZfKac/JvXucP/GEeGoICBBPC3Pm5H7N+fOcv/ce588+y/mMGeIGkzPDRyazzxvXaDhv1kzUH1Svbi9RoNeLXPYvvpDWsl+woOC/UXZsNs737OF88WLOT58umjEJogRQOhz+449LO/yAAM7378/fWN7g2Wedr6aVSs6Dgji/cMG1sa5dE6X7ed2A/vxTOOGM8IhOx3nlyo5PSpGRnE+fLjRf6tYVWvwpKeIG4SwE5EwTv3r1vO3ftk2knrZsKdIXHz507XsTBFFK0jIHDZLetGNMtMbzdXr2FFLFUpjNQsVx7Fjn1ycnA2+/LVIlIyKAH3+0V8zMic0GvPSSyOrJkGBISxOZP0aj/bkxMUL6+ORJ4OxZoaKp1wPbtzvPCnKWwXT1KnDpknO7Jk0SCqObNonUz88/F9k7OVVBCYIoFMXb4Y8cKfRfMpymQiEyXxYuzF9Tbm8xZAhQs6ZIUZTCZgN27HD+2RNPiN63CQkiVXPRIiFX7CwN9MIF6RuCzeborA0GcQPJSfXquWcGScE50KWLdJ1AYiLw3//a30QMBuDGDdd77BIE4RLF2+HrdMCBA8Ds2WK1/847oqFFv37etsw1rFaRxvjWW85TMQMCpI/v3CkceHYRMLNZFCOtXWt/7tmz4kng7belBcyckZLiWDD15psFu5nGxYnVe04OHZJOuzQYRAonILp9DR4MdOwofi9a+RNEgSj+efhqNfDii+JVnJg1SzQkUShEeCUoSKy+s6/OdTpxE7PZRKHV7NliJfzCC+JmIbWST0kBjh7NKoDauFH82WQS8+RHgbJpU7ECj4sTTyJqtXiiWrUKePllURAltWqXQiYD7t1zPF6+vHQoiDGgYkWhBjpmjLgBcC5uGj/9lFW3QBCE67gS6PfUy+e1dIqKXbuk1Sb9/e1VHIcPF+mcw4fbb6pqtWIjVEqfRq/PyoqxWISgmlSOfm5Nx+VyYd+TT4pMnYwah5kzs76D1SpULnMKqSmV0oJsGg3n8fGOv4XNJjaFs4u+ZWwm79wpvUGs0Yjcf4IgOOeub9oW/xW+N7DZhIa7UgnUrp3/toXff++48ZkRR1+zRlQKN2ggwjkvvQT89pv9uQaDWF1rtfbt+mQy8VSQsbo/e9ZxMxYQblOtdr752rKlmHv3bhECyggDjR0LVKsmWibKZMBXX4nvPn26eFIxm4W9e/YInfyMymK9XrQQLF/ecS7GgK1bgWeeEfIJCoX4Pv/7n/gdpGQHjMYsbSGCIFyGHH5++ftvsV+QnCwcZ5UqQs+lfn3Xx7hzR/q4QiGcY/v2YuyWLUV4RorUVOEkHz4UUgcA0KaN2LDO2MT297dviJIdjca5wz95Ujj5nCGjtDTg66+FwweEs/7qK+DTT4UuTJUqYs6UFCEYt2qVyCB6+23H/rbZqVZNzHn6tAghNWsmblwnTzq3v0IF5+MRBCGNK48Bnnr5fEjn1i3pyt7gYCEn4IzYWCHLkJAg3k+cKK0pr9WKBh2cc/7337lXwKpUnE+YIM41Gp3n37do4Ri+0euFJK8zWQaFwrlOflhY0f2eeWGzcd6woXS4Z8eO/I/1yy+irqBSJSHLfPWqW8wmCE+DUpGH72l+/dVxxcm5CDFs2iTeX7wo2uU99ZQIObRtKwTR+vQBKlUSq+E33wRCQuzTMXU6kY+eIdB28qTzvHZAPA288or4s1rtPLVzzRqgVi0xrr+/OO+VV4DJk8UThBQWi7ScsUwmUkE9BWPA5s1AnTriqSUgQNj/+edA5875G2vsWPG7nzolRNyWLBFPEs6etgiiJOLKXcFTL59e4ZtMnPfpI73qVSo5r1FDNPBWqbJWxzlXphmr6yVLROOSiRNFtXDv3mKDMjvbtjlf4ev1QqzMVWw2oZmzejXn169nHZ82zfkTRLVq9humcrmoYI6NLZKfM1/YbJwfOSJ+k5zN0V3h3j3pJyqVSmgCEUQxB6VCWsFTREVxXqaM8xBIfl9NmuQ9p9XKeb169qEVmUxk8GSEhgrLokXSNyWA8wEDxE3liSfEzezll3Nv4+jL7N4tfjep79mypbetI4hC46rDp03bvEhOFpujyclFN6ZUPnpOZDKR7TJ6tNgU5hzo0EFU1pYpUzR2dOoksmByFmNptULDv2tX8SruVKkiXbMgk4lwF0GUEiiGnxfr1glnmxPGRJu6/MoMKBS5Z6xkp3x50SkqPV3sE+zcWbQOqnJl0XVKp8s6pteLG8szz9ifm54OzJghGne0awcsXiwd5/dFatUCWrd2rBDWaID33/eOTQThBWiFnxcPH0pvnnIuHF9UlOurf7VabJz+97/5s0Emy1+FbH4YPx548kmhW5OSInL4+/e37xpmtYqngWPHslI5jx0TbfoWL3aPXUXNmjXA8OFiE1gmE5XNc+YUD5E9gigiyOHnxVNPSR/X60XWx7590p/7+4tCpCefFNklV64ILZgxY0SGjq+wZ4/IWMnI3iJ+pp4AAA5NSURBVOnQwbGQbONG0Qg9e95+aiqwerXIfmnUyLM2F4TAQOH0Hz4UrypV3HcTJQgfhRx+XtStC7z6KvDzz1lKk3q9SAvs1g3YtUukXN67JxyIVivSK4ODhSOsVs279ufG++8LXZqMitilS8UqeOZM+/P++su5YFlUlL3D51wofC5aJEI+Q4eK3ym/1cjuIjBQvAiiNOLKzq6zF4ABAGIA2ABE5vjsEwCxAM4BeNqV8Xw2S8dm43zzZs779+e8Vy/Oly8XOjXZPz9xQmTzjBsnWixGRHD+44/25/kSJ09KZx1ptZwfPWp/7tdf23e/ynj5+3O+cqX9uW++aV+cptcLLSCCINwGPJGWCaA+gLoAdmd3+AAaADgOQA2gBoCLAOR5jeezDt8VjEaRh58931un43zQIG9bJs2kSdLVtHI5519+aX/ujRvSImblytlX+J44IX2eTsf5v/969vsRRCnCVYdfqCAm5/wM5/ycxEd9ACzjnKdzzi8/Wuk/Xpi5fJ4VK4SeTHaxsrQ0keVz+rT37HKGTiedYZTRRCY7lSsLbfry5cXehFYrztFogOeey9K5//NPsW+RE6MR2LKl6GznHLh+XVTMEgThMu7ataoM4Hq29zceHSu57Nwp3U1KJhNNWnyN556TPs4YMHCg4/HOnYUMwcyZwuGmpQE3bwpJiU6dRMaOv7+0uqVSWXTa9UePCqG6OnVEumWTJsA5qTUHQRA5ydPhM8Z2MMZOSbz65HaZxDGJZHaAMTaCMXaYMXY4Pj7eVbt9j+rVpTs3yWRCQ8fXCA0VqZgKRVbap1Ip1DYrO7k3y+UilTGn5HJamlDEdHYTkcmyJJsLw4MHItPp3Dlhg9Eosofat5eWgSYIwo48HT7n/CnOeSOJ17pcLrsBoGq291UASD5/c87ncs4jOeeRwcHB+bPel3j1VccQiUwmMkK6dPGOTblhMgHffSecvM0mXkqlyDrKjSNHpI9fuCDST1euzBI6CwgQf166tGhuer//7hgy4o/E69bl9teRIAjAfSGd9QAGMcbUjLEaAGoD+NdNc/kGVaqIop4qVbLi202aiLTF7EVMvsLq1cD581kpmYBYqS9eDMTGOr9OqokJIBy7Ugn06CFaIi5ZIsa6exfo27dobL5+XVrD32gUTc8JgsiVQjl8xlhfxtgNAK0BbGKMbQUAznkMgBUATgP4E8CbnPNctH5LCB06iI3bY8fEijc6GqhRw71zJiYCb7whHHGFCqKwy1nOfHaJiG3bpM+Ty0UxljM++MBxUxcQq/mEBPFnnU40SendO6sZS1HQqlWWfHR2VCrnUs8EQWRS2CydtZzzKpxzNec8hHP+dLbPvuKc1+Kc1+WcF2GKho/DmGh7WKWK++cym4VGzIIFwP37ou3h7Nmiuje7zs3OnaI4Si4HypYFvvxShFhyassA4pzcKoHffVds0uYkLg4YMqTQXylXnnlGbNZm1/7XaoW+T9u27p2bIEoAVFtenFm3TmTKZFeCTE8Xm5oZbQ///Rfo2ROIiREr/MRE4JtvhIPOuefAmHCgue05MCbmzInFInrguqIEWlDkchEiGztWZOjUrg385z8i5dNXKnkJwochh1+cOXJEOixjNIqwEiC6Q2WP0wMiDr50qYizlysnUib1euFA//5bOrUyOxmhm5zI5UKnxp3o9UJ8LjZW7EGMHSudHUUQhAPk8Isz4eHSMXKtFqhZU/w5Jkb6WrkcaNBA5Nbv3g0cOgScPSu0g/Kie3fpoi0/P5GeShCET0IOvzgzcKBw7tlVH+VysYHaq5d437ixdLjDahX7DAqFkAiuX9/1sMiECWIvICOWLpOJjdp583wzI4kgCADk8Is3fn7A/v1i41apFK8OHYRkc8aG7H/+I24K2dHpgHfeKXgGTcWK4snh449FT4ChQ4F//hFZOQRB+CyMS3Vz8hKRkZH88OHD3jajeJKSIlboUk48Kkqkax4/LmL2H3wgpJFJD54gSgSMsWjOeWRe55EefklBKj89gw4dRE0AQRClGlriEQRBlBJohU/4Flar0PO5fVvsTYSHe9sigigxkMMnpElNBVatEr14IyNFm0J3Z+BcuSKqhBMSRJGYxQIMGiQqiWm/gSAKDTl8wpFz50T2jdEoNoP9/MRKe8+e3PcKCku/fkIgLbssxIoVQv74lVfcNy9BlBJo2UQ4MnSo0ObJqOJNSQHOnBEaPO7i6lUxR3ZnD4iq4B9/dN+8BFGKIIdP2HP/vmgqkjNdNz1dyB27C4PBecjImfonQRD5ghw+YU9u1bbuFCirU0e6DaJGI91ykSCIfEMOn7CnbFnRuCXnJqlGA7z4ovvmlcmEmJtOl1UlrNcLbZ4PPnDfvARRiqBNW8KRJUvEpm1amnjpdEC9esBnn7l33s6dhWTDvHkiY6dLF5Glk13/niCIAkPSCsUZo1FIBS9YIGLs3bsDU6cCVavmeWmeGAzAmjViMzUyEnjqKUqNJAgfxVVpBXL4xZmnnxY6OUajeC+XC62cc+eAoCDv2kYQhMdw1eHTkq24cvw4sHdvlrMHRJVqSgrwyy9eM4sgCN+FHH5x5cQJ6RBLWhpw4IDn7SEIwuchh19cqV3bMVceEBucjRt73h6CIHwecvjFlZYtRTvCjBTGDFQq4LXXvGMTQRA+DTn84gpjwI4dQJ8+otOVXA48/riI64eEeNs6giB8EMrDL86UKSPExcxmsWFL+eoEQeQCOfySQEY/W4IgiFygkI4vkpICTJggJInr1QOmTAFMJm9bRRBEMYdW+L6GxSL038+ezcqx/89/RLz+zz/dK2BGEESJhlb4vsaGDUBsrH1BlcEA/PMPcPCg9+wiCKLYUyiHzxgbwBiLYYzZGGOR2Y53YYxFM8ZOPvp3p8KbWkrYt09a/91sJodPEEShKGxI5xSAfgDm5Dh+D0AvzvktxlgjAFsBVC7kXKWDqlWFOmVamv1xtRqoTD8hQRAFp1ArfM75Gc75OYnjRznntx69jQGgYYypCzNXqWHoUECR4z7MGKDVAr16eccmgiBKBJ6I4fcHcJRznu6BuYo/ZcsCO3cK6QStVuTWN2kiGoir6Z5JEETByTOkwxjbASBU4qNxnPN1eVzbEMAkAF1zOWcEgBEAEBYWlpc5pYPmzYXE8fXrYrVfqZK3LSIIogSQp8PnnD9VkIEZY1UArAXwEuf8Yi7jzwUwFxB6+AWZq0TCGEA3QIIgihC3hHQYY0EANgH4hHP+jzvmIAiCIPJHYdMy+zLGbgBoDWATY2zro4/eAhAOYDxj7NijV4VC2koQBEEUgkKlZXLO10KEbXIe/xLAl4UZmyAIgihaqNKWIAiilEAOnyAIopTAuFSbPC/BGEsG4FDI5eOUh6gsLm4UR7vJZs9RHO0uzTZX45wH53WSr6llnuOcR+Z9mu/AGDtc3GwGiqfdZLPnKI52k815QyEdgiCIUgI5fIIgiFKCrzn8ud42oAAUR5uB4mk32ew5iqPdZHMe+NSmLUEQBOE+fG2FTxAEQbgJn3T4jLG3GWPnHnXTmuxte1yFMfYBY4wzxsp725a8YIxNYYydZYydYIytfaR/5JMwxro9+vsQyxj72Nv2uAJjrCpjbBdj7Myjv8fvetsmV2GMyRljRxljG71tiyswxoIYY6se/X0+wxhr7W2bXIExNubR341TjLHfGWMad8/pcw6fMdYRQB8AEZzzhgCmetkkl2CMVQXQBcA1b9viItsBNOKcRwA4D+D/27ufEKvKOIzj3wcnKwvBTRsnGAecUpJwSJGkCMdFf2TWBoboKjE1CDKN9q4soXY6bmYgchqkRX8MgnZaNCaTziYMdPyDszEFF8Mwj4v3CAOO5x7Te99zu7/Pas5h4DwwZ5573vfcc94DmfMsSNIi4CvgLWA18K6k1XlTVTILfGR7FbAB2N0muQH2AZO5QzyEI8CPtl8EXqYNsktaDuwFXrH9ErAI2Nrs49au8IFdwKF7C6bYvpE5T1WfAx8DbXFTxPYp27PF5mmgO2eeEuuBv21ftD0DfE26IKg129dsjxc/3yaVUO3XqCxea/4OcDR3liokLQVeB44B2J6xfTNvqsq6gKcldQFLgKsNfv+R1bHw+4DXJJ2R9KukdbkDNSJpELhi+1zuLP/RTuCH3CEeYDlwed72FG1QnPNJ6gHWAu2wCv0XpAuXudxBKuoFpoHjxTTUUUnP5A7ViO0rpNmLS8A14F/bp5p93CxP2patokXKtIw0DF4HfCOp15m/TtQg80FKVvXKpcpqZZI+JU0/jLQy20PQAvvaYhQFIOlZ4FvgQ9u3cucpI2kLcMP2H5LeyJ2noi6gH9hj+4ykI8AnwGd5Y5WTtIw0Ul0B3AROSNpme7iZx81S+GWraEnaBYwVBf+bpDnS+yamW5VvIQ/KLGkN6Y92ThKkqZFxSettX29hxPs0Wq1M0nZgCzCQ+wO1xBTw/Lztblow9H0cJD1BKvsR22O581SwERiU9DbwFLBU0rDtbZlzlZkCpmzfGz2Nkgq/7jYD/9ieBpA0BrwKNLXw6zilcxLYBCCpD1hMjV+IZHvC9nO2e2z3kE7A/txl34ikN4H9wKDtO7nzlPgdWClphaTFpBtb32XO1JDSp/8xYNL24dx5qrB9wHZ3cR5vBX6pedlT/J9dlvRCsWsAuJAxUlWXgA2SlhTnygAtuNlct5enAQwBQ5L+AmaA7TW++mxnXwJPAj8XI5PTtt/PG+l+tmclfQD8RPomw5Dt85ljVbEReA+YkPRnse+g7e8zZvq/2gOMFBcEF4EdmfM0VEw/jQLjpCnVs7Tgqdt40jaEEDpEHad0QgghNEEUfgghdIgo/BBC6BBR+CGE0CGi8EMIoUNE4YcQQoeIwg8hhA4RhR9CCB3iLtxoppTS+NmWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "# 导入sklearn模拟二分类数据生成模块\n",
    "from sklearn.datasets.samples_generator import make_blobs\n",
    "# 生成模拟二分类数据集\n",
    "X, y =  make_blobs(n_samples=150, n_features=2, centers=2,\n",
    "  cluster_std=1.2, random_state=40)\n",
    "# 将标签转换为1/-1\n",
    "y_ = y.copy()\n",
    "y_[y_==0] = -1\n",
    "y_ = y_.astype(float)\n",
    "# 训练/测试数据集划分\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y_,\n",
    " test_size=0.3, random_state=43)\n",
    "# 设置颜色参数\n",
    "colors = {0:'r', 1:'g'}\n",
    "# 绘制二分类数据集的散点图\n",
    "plt.scatter(X[:,0], X[:,1], marker='o', c=pd.Series(y).map(colors))\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 定义决策树桩类\n",
    "### 作为Adaboost弱分类器\n",
    "class DecisionStump():\n",
    "    def __init__(self):\n",
    "        # 基于划分阈值决定样本分类为1还是-1\n",
    "        self.label = 1\n",
    "        # 特征索引\n",
    "        self.feature_index = None\n",
    "        # 特征划分阈值\n",
    "        self.threshold = None\n",
    "        # 指示分类准确率的值\n",
    "        self.alpha = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 定义AdaBoost算法类\n",
    "class Adaboost:\n",
    "    # 弱分类器个数\n",
    "    def __init__(self, n_estimators=5):\n",
    "        self.n_estimators = n_estimators\n",
    "        \n",
    "    # Adaboost拟合算法\n",
    "    def fit(self, X, y):\n",
    "        m, n = X.shape\n",
    "        # (1) 初始化权重分布为均匀分布 1/N\n",
    "        w = np.full(m, (1/m))\n",
    "        # 处初始化基分类器列表\n",
    "        self.estimators = []\n",
    "        # (2) for m in (1,2,...,M)\n",
    "        for _ in range(self.n_estimators):\n",
    "            # (2.a) 训练一个弱分类器：决策树桩\n",
    "            estimator = DecisionStump()\n",
    "            # 设定一个最小化误差\n",
    "            min_error = float('inf')\n",
    "            # 遍历数据集特征，根据最小分类误差率选择最优划分特征\n",
    "            for i in range(n):\n",
    "                # 获取特征值\n",
    "                values = np.expand_dims(X[:, i], axis=1)\n",
    "                # 特征取值去重\n",
    "                unique_values = np.unique(values)\n",
    "                # 尝试将每一个特征值作为分类阈值\n",
    "                for threshold in unique_values:\n",
    "                    p = 1\n",
    "                    # 初始化所有预测值为1\n",
    "                    pred = np.ones(np.shape(y))\n",
    "                    # 小于分类阈值的预测值为-1\n",
    "                    pred[X[:, i] < threshold] = -1\n",
    "                    # 2.b 计算误差率\n",
    "                    error = sum(w[y != pred])\n",
    "                    \n",
    "                    # 如果分类误差大于0.5，则进行正负预测翻转\n",
    "                    # 例如 error = 0.6 => (1 - error) = 0.4\n",
    "                    if error > 0.5:\n",
    "                        error = 1 - error\n",
    "                        p = -1\n",
    "\n",
    "                    # 一旦获得最小误差则保存相关参数配置\n",
    "                    if error < min_error:\n",
    "                        estimator.label = p\n",
    "                        estimator.threshold = threshold\n",
    "                        estimator.feature_index = i\n",
    "                        min_error = error\n",
    "                        \n",
    "            # 2.c 计算基分类器的权重\n",
    "            estimator.alpha = 0.5 * np.log((1.0 - min_error) / (min_error + 1e-9))\n",
    "            # 初始化所有预测值为1\n",
    "            preds = np.ones(np.shape(y))\n",
    "            # 获取所有小于阈值的负类索引\n",
    "            negative_idx = (estimator.label * X[:, estimator.feature_index] < estimator.label * estimator.threshold)\n",
    "            # 将负类设为 '-1'\n",
    "            preds[negative_idx] = -1\n",
    "            # 2.d 更新样本权重\n",
    "            w *= np.exp(-estimator.alpha * y * preds)\n",
    "            w /= np.sum(w)\n",
    "\n",
    "            # 保存该弱分类器\n",
    "            self.estimators.append(estimator)\n",
    "    \n",
    "    # 定义预测函数\n",
    "    def predict(self, X):\n",
    "        m = len(X)\n",
    "        y_pred = np.zeros((m, 1))\n",
    "        # 计算每个弱分类器的预测值\n",
    "        for estimator in self.estimators:\n",
    "            # 初始化所有预测值为1\n",
    "            predictions = np.ones(np.shape(y_pred))\n",
    "            # 获取所有小于阈值的负类索引\n",
    "            negative_idx = (estimator.label * X[:, estimator.feature_index] < estimator.label * estimator.threshold)\n",
    "            # 将负类设为 '-1'\n",
    "            predictions[negative_idx] = -1\n",
    "            # 2.e 对每个弱分类器的预测结果进行加权\n",
    "            y_pred += estimator.alpha * predictions\n",
    "\n",
    "        # 返回最终预测结果\n",
    "        y_pred = np.sign(y_pred).flatten()\n",
    "        return y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy of AdaBoost by numpy: 0.9777777777777777\n"
     ]
    }
   ],
   "source": [
    "# 导入sklearn准确率计算函数\n",
    "from sklearn.metrics import accuracy_score\n",
    "# 创建Adaboost模型实例\n",
    "clf = Adaboost(n_estimators=5)\n",
    "# 模型拟合\n",
    "clf.fit(X_train, y_train)\n",
    "# 模型预测\n",
    "y_pred = clf.predict(X_test)\n",
    "# 计算模型预测准确率\n",
    "accuracy = accuracy_score(y_test, y_pred)\n",
    "print(\"Accuracy of AdaBoost by numpy:\", accuracy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy of AdaBoost by sklearn: 0.9777777777777777\n"
     ]
    }
   ],
   "source": [
    "# 导入sklearn adaboost分类器\n",
    "from sklearn.ensemble import AdaBoostClassifier\n",
    "# 创建Adaboost模型实例\n",
    "clf_ = AdaBoostClassifier(n_estimators=5, random_state=0)\n",
    "# 模型拟合\n",
    "clf_.fit(X_train, y_train)\n",
    "# 模型预测\n",
    "y_pred_ = clf_.predict(X_test)\n",
    "# 计算模型预测准确率\n",
    "accuracy = accuracy_score(y_test, y_pred_)\n",
    "print(\"Accuracy of AdaBoost by sklearn:\", accuracy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
